
COMPILE_OPT IDL2, LOGICAL_PREDICATE

COMMON OMNI_CONFIG, conf, mw, local, dpdfs, ancil, fmt, conffile
COMMON HRDS_BLOCK, hrds, hdist
COMMON PARALLAX_BLOCK, bessel, adist
COMMON VEL_BLOCK, v, v_std


;; Read in galactic-params & dpdf-params config files
conf = omni_load_conf()
IF ~exist(mw) THEN mw = omni_read_conffile('./conffiles/galactic_params.conf')
IF ~exist(ancil) THEN ancil = omni_read_conffile('./conffiles/ancillary.conf')
IF ~exist(dpdfs) THEN $
   dpdfs = omni_read_conffile('./conffiles/dpdf_params.conf')

s = omni_read_cat(conf.cat,ncat)
restore,'./local/'+conf.survey+'_pvec.sav',/ver
restore,'./local/'+conf.survey+'_velocities.sav',/ver
v_std = findgen(conf.nvbin)*conf.deltav + conf.vstart
d = findgen(dpdfs.nbins) * dpdfs.binsize + dpdfs.binstart

;; Select the subset of the BGPS that reasonably overlaps with the
;;   HRDS region for computation -- speeds things up nicely.
lind = where(s.glon GE 13. AND s.glon LE 70. AND $
             (constrain.kdist OR constrain.grsmatch OR $
              constrain.parallax OR constrain.knownd), nl)
print,'Number in culled sample: ',nl

s         = s[lind]
pvec      = pvec[lind]
constrain = constrain[lind]
v         = v[lind]
ncat      = nl

ndr = fix(floor(4*5*1.5))
ndv = fix(floor(5*5*1.5))


dr = (dindgen(ndr)*1.d + 1.d)*1.d
dv = dindgen(ndv)*0.2d + 2.d


goto,plot_only


blank = lonarr(ndr,ndv)

nkd = blank
ngr = blank
nem = blank
nh2 = blank
nkn = blank
npx = blank
nhr = blank
npo = blank

ndi = blank

nma = blank
fma = double(blank)



start_t = systime(1)

dcomp = replicate({hrds:-1000.d,$
                   post:-1000.d},ncat)


nloop = 0

FOR ii=0,ndr-1 DO BEGIN
   FOR jj=0,ndv-1 DO BEGIN
            
      message,'Now starting loop '+$
              string(++nloop,ndr*ndv,format="(I0,' of ',I0)")+'...',/inf
      
      ;; Clear array
      dcomp.hrds = -1000.
      dcomp.post = -1000.
      
      ;; Set parameter values for this run
      ancil.ppv_dis = dr[ii]
      ancil.ppv_dv  = dv[jj]
      undefine,adist,hdist
      
      ;; Run!
      ;; distance_omnibus
      FOR kk=0,ncat-1 DO BEGIN
         
         ;;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
         ;; Set flag for rejection of kinematic information based on
         ;;   location in the L-V diagram
         ;; Choose VLSR from dense gas or GRS 13CO
         CASE 1 OF
            v[kk].vlsr GE -500.: vlsr = v[kk].vlsr                   
            where([1,2,5,6] EQ v[kk].grs.flag) NE -1: vlsr = v[kk].grs.vlsr 
            ELSE: vlsr = -1000.
         ENDCASE
         reject_kd = OMNI_KINEMATIC_AVOIDANCE(v[kk].l,vlsr)
         ;;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
         
         p_hrds = prob_hrds(s[kk], CONSTRAIN=c_hrds)
         constrain[kk].hrds     = c_hrds
         
         ;;=================================================================
         ;; Check which DPDFs to run
         ntags = n_tags(dpdfs.dpdf)
         dvals = intarr(ntags)
         FOR tag=0,ntags-1 DO dvals[tag] = dpdfs.dpdf.(tag)
         dpdfi = where(dvals, nrun)
         IF ~nrun THEN BEGIN
            message,'Error: No DPDFs selected to run in dpdf_params.conf!  '+$
                    'Exiting.',/cont
            RETURN
         ENDIF
         ;; DTAGS contains only the DPDF tags in the order specified
         dtags = (tag_names(dpdfs.dpdf))[dpdfi[sort(dvals[dpdfi])]]
         has_kin = fix(total(strmatch(dtags,'KDIST'))) || $
                   fix(total(strmatch(dtags,'GRSMATCH')))
         ;;=================================================================
         ;;=================================================================
         ;; Evaluation of the posterior DPDF
         ;;   Criteria are laid out in Ellsworth-Bowers et al. (2013, ApJ,
         ;;   770, 39).
         
         dtan = mw.r0 * cos(s[kk].glon*!dtor) / cos(s[kk].glat*!dtor)
         pvec[kk].stat.dtan = dtan
         
         kine = constrain[kk].kdist ? $
                pvec[kk].kdist : pvec[kk].grsmatch
         
         ;;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
         ;;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
         ;; Repeat 2x for kine + EMAF + H2  && kine + HRDS
         FOR hh=0, 1 DO BEGIN
            
            constrain[kk].post = 0b ; Set to zero for each case
            
            priors = hh ? p_hrds : (pvec[kk].h2*pvec[kk].emaf*pvec[kk].knownd)
            pvec[kk].post = kine * priors
            pvec[kk].post /= total(pvec[kk].post)
            
            ;; Compute D_ML and associated error bars
            pvec[kk].stat.dml = OMNI_COMPUTE_EBARS( pvec[kk].post, d, $
                                                    PVAL=0.6827d )
            
            ;; print,pvec[kk].stat.dml
            ;; Compute FW68
            pvec[kk].stat.fw68 = pvec[kk].stat.dml[1] + pvec[kk].stat.dml[2]
            
            ;; If INNER GALAXY, find regions on either side of DTAN
            IF (s[kk].glon GT -90. && s[kk].glon LT 90.) || $
               (s[kk].glon GT 270.) THEN BEGIN
               
               nind = where(d LE dtan, nnear)
               find = where(d GT dtan, nfar)
               ;; Compute P_ML
               pvec[kk].stat.pml = total(pvec[kk].post[nind]) > $
                                   total(pvec[kk].post[find])
               ;; Make N/F/T determination -- Tangent if <= 1 kpc of dtan
               pvec[kk].stat.kdar = pvec[kk].stat.dml[0] LE dtan ? 'N' : 'F'
               IF abs(pvec[kk].stat.dml[0] - dtan) LE 1.d3 THEN $
                  pvec[kk].stat.kdar = 'T'
               
            ENDIF ELSE BEGIN
               ;; OUTER GALAXY sources, P_ML = 1., KDAR = 'O'
               pvec[kk].stat.pml = 1.d
               pvec[kk].stat.kdar = 'O'
            ENDELSE
            
            
            ;; Compute moment distances (dbar)
            pvec[kk].stat.dbar[0] = total(pvec[kk].post * d)
            pvec[kk].stat.dbar[1] = sqrt( total(pvec[kk].post * d * d) - $
                                          pvec[kk].stat.dbar[0] * $
                                          pvec[kk].stat.dbar[0] )
            
            
            ;; Fill in DUSE elements
            pvec[kk].stat.duse = (pvec[kk].stat.kdar EQ 'T') ? $
                                 [pvec[kk].stat.dbar[0],$
                                  pvec[kk].stat.dbar[1],$
                                  pvec[kk].stat.dbar[1]] : $
                                 pvec[kk].stat.dml
            ;;=================================================================
            ;;=================================================================
            
            
            
            ;;;;;;;;;;;;;;;;
            ;; Check to see if survey kas kinematic information.  If not,
            ;;   then skip this step.  This check is to allow the code to not
            ;;   puke on SURVEYS that have no kinematic information.
            IF has_kin THEN BEGIN
               
               ;;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
               ;; Set POSTERIOR CONSTRAINT and KDAR flags based on statistics
               ;;   values and other relevant information.
               
               ;; To have a "well-constrained" distance estimate, a SURVEY
               ;;   source must have a constraint from at least ONE of the
               ;;   following:   [KDIST, GRSMATCH, KNOWND, PARALLAX]
               constrain_arg = constrain[kk].kdist
               IF conf.usegrs THEN $
                  constrain_arg = constrain_arg || constrain[kk].grsmatch 
               IF dpdfs.dpdf.knownd NE 0 THEN $
                  constrain_arg = constrain_arg || constrain[kk].knownd
               IF dpdfs.dpdf.parallax NE 0 THEN $
                  constrain_arg = constrain_arg || constrain[kk].parallax
               
               ;; Set CONSTRAIN[kk].POST=1b IF above criteria are met, AND the
               ;;   FW68 is less than 2.3 kpc *OR* D_ML is w/in 1 kpc of DTAN
               IF constrain_arg && $
                  (pvec[kk].stat.fw68 LE 2.3d3 || $
                   abs(pvec[kk].stat.dml[0] - pvec[kk].stat.dtan) LE 1.d3) $
               THEN $
                  constrain[kk].post = 1b ELSE $
                     pvec[kk].stat.kdar = 'U' ; If CONSTRAIN[kk].post = 0
               
               ;;**************************************************************
               ;; Check for whether this source is in the Galactic bar --
               ;;   don't serve it any more alcohol -- also, reject
               ;;   sources near l=90 and l=180, unless **KNOWND OR PARALLAX**
               ;;   are set.  Essentially, looking at the KINEMATIC-ONLY
               ;;   sources that passed the above test, and setting KDAR='X'
               ;;   if they are in a rejected region.  However, this code
               ;;   should NEVER be run, since we already set kinematic
               ;;   constraints to zero above...
               
               has_kid = 0b     ; Kinematic-Independent Distance Measurement
               IF dpdfs.dpdf.knownd NE 0 THEN $
                  has_kid = has_kid || constrain[kk].knownd
               IF dpdfs.dpdf.parallax  NE 0 THEN $
                  has_kid = has_kid || constrain[kk].parallax
               
               IF reject_kd && ~has_kid THEN BEGIN
                  ;; Set the CONSTRAIN.POST and KDAR correctly if adjusted above
                  constrain[kk].post = 0b
                  pvec[kk].stat.kdar = 'X' ; X for eXcluded based on L-V
               ENDIF
               
               ;;************************************************************
               ;;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
            ENDIF
            ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
            
            
            IF constrain[kk].post THEN BEGIN
               CASE hh OF
                  0: BEGIN
                     dcomp[kk].post = pvec[kk].stat.duse[0]
                     IF dcomp[kk].post LE 100 THEN BEGIN
                        cc = constrain[kk]
                        print,cc.kdist,cc.grsmatch,cc.h2,cc.emaf,cc.hrds
                        
                        cgPlot,d,pvec[kk].emaf,color='blue'
                        cgOplot,d,kine,color='forest green'
                        cgOplot,d,pvec[kk].post,color='black'
                        wait,0.5
                     ENDIF
                  END
                  1: IF constrain[kk].hrds THEN $
                     dcomp[kk].hrds = pvec[kk].stat.duse[0]
               ENDCASE
            ENDIF
         ENDFOR                 ; End of HRDS / POST mini-loop
         ;;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
         ;;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
         
         ;; line = 0b
         ;; plin = 0b
         ;; IF constrain[kk].hrds THEN BEGIN
         ;;    print,'CH: ',constrain[kk].hrds,dcomp[kk].hrds
         ;;    line = 1b
         ;;    plin++
         ;; ENDIF
         ;; IF dcomp[kk].hrds GT 0. THEN BEGIN
         ;;    print,'DHRDS: ',constrain[kk].hrds,dcomp[kk].hrds
         ;;    line = 1b
         ;;    plin++
         ;; ENDIF
         ;; IF line THEN BEGIN
         ;;    IF plin NE 2 THEN BEGIN
         ;;       cgPlot,d,pvec[kk].post
         ;;       cgOplot,d,pvec[kk].hrds,color='orchid'
         ;;       print,pvec[kk].stat.kdar
         ;;       wait,0.5
         ;;    ENDIF
         ;;    print,'-------------------------------'
         ;; ENDIF
         
      ENDFOR                    ; End of loop over sources
      
      
      ;; Load up arrays
      nkd[ii,jj] = long(total(constrain.kdist))
      ngr[ii,jj] = long(total(constrain.grsmatch))
      nem[ii,jj] = long(total(constrain.emaf))
      nh2[ii,jj] = long(total(constrain.h2))
      nkn[ii,jj] = long(total(constrain.knownd))
      npx[ii,jj] = long(total(constrain.parallax))
      nhr[ii,jj] = long(total(constrain.hrds))
      npo[ii,jj] = long(total(constrain.post))
      
      ;; Compare distances
      ind = where(dcomp.hrds GT 0. AND dcomp.post GT 0., nind)
      dd = (dcomp[ind].hrds - dcomp[ind].post) / 1.d3
      
      ;; Little plot for user's amusement
      plothist,dd,bin=0.1,xarr,yarr,xtit=cgSymbol('Delta')+'d  [kpc]'
      yf = mpfitpeak(xarr,yarr,A,nterms=3)
      cgOplot,xarr,yf,color='red',thick=3
      print,A
      
      matchi = where(abs(dd) LE 1., nmatch) ; <= 1 kpc == good
      print,'N HRDS / N dist / N match: ',nhr[ii,jj],nind,nmatch
      
      ;; Load up arrays
      ndi[ii,jj] = nind
      nma[ii,jj] = nmatch
      fma[ii,jj] = double(nmatch) / nind
      
      al_legend,/top,/right,box=0,$
                ['Loop '+string(nloop,ndr*ndv,format="(I0,'/',I0)"),$
                 'R = '+string(ancil.ppv_dis,format="(F0.1,' pc')"),$
                 'V = '+string(ancil.ppv_dv,format="(F0.1,' km/s')"),$
                 'Dis / Match / f: '+string(ndi[ii,jj],nma[ii,jj],fma[ii,jj],$
                                            format="(I0,' / ',I0,' / ',F0.3)")]
      
      
      ;; Figure time elapsed and ETD...
      elapsed = systime(1) - start_t
      done = double(nloop) / double(ndr*ndv)
      print,'ETD: ',systime(0,start_t + elapsed/done),$
            '  ['+string(nloop,ndr*ndv,format="(I0,'/',I0,']')")
      
      
      ;; Save at each step along the way!
      save,nkd,ngr,nem,nh2,nkn,npx,nhr,npo,ndi,nma,fma,$
           filename='./grs_paper/assoc_volume.sav'
      
   ENDFOR
ENDFOR




;;()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()
;;()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()()
plot_only:

restore,'./grs_paper/assoc_volume.sav',/ver


myps,'./grs_paper/assoc_volume.eps',xsize=10,ysize=10

multiplot_xm,mpcharsize=1.0,[2,2],gap=0.04,/doxaxis,/doyaxis

contcol = 'blk2'
valcol = 'black'

;;================================
;; N PARALLAX

print,m4_stat(npx)
cgLoadct,13
pr = set_plot_range(npx)
print,'N PARALLAX',pr
plotimage,npx,xst=4,yst=4,range=pr,title='Maser Parallax',min_dpi=100

ddr = mean(dr[1:*] - dr[0:*])
ddv = mean(dv[1:*] - dv[0:*])

cgAxis,xaxis=0,xr=[min(dr),max(dr)]+[-1,1]*ddr/2.,/xst,charsize=1.0,/save,$
       xtit='Physical Radius  [pc]'
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
cgAxis,yaxis=0,yr=[min(dv),max(dv)]+[-1,1]*ddv/2.,/yst,charsize=1.0,/save,$
       ytit='FWHM Velocity Dispersion  [km s!u-1!n]'
cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

cgContour,npx,dr,dv,levels=findgen(21)*50.,/over,color=contcol

cgPlots,psym=14,color='brown',symsize=2,  16.5, 5.1
cgPlots,psym=4,color='black',symsize=2,   16.5, 5.1
cgPlots,psym=14,color='bisque',symsize=2,  7.4, 3.1
cgPlots,psym=4,color='black',symsize=2,    7.4, 3.1
cgPlots,psym=14,color='blu3',symsize=2,   14.7, 4.75
cgPlots,psym=4,color='black',symsize=2,   14.7, 4.75
cgPlots,psym=14,color='pink',symsize=2,   10.5, 3.6
cgPlots,psym=4,color='black',symsize=2,   10.5, 3.6
cgPlots,psym=14,color='grn4',symsize=2,    9.0, 4.1
cgPlots,psym=4,color='black',symsize=2,    9.0, 4.1

multiplot,/doxaxis,/doyaxis

;;================================
;; N HRDS

print,m4_stat(nhr)
cgLoadct,13
pr = set_plot_range(nhr)
print,'N HRDS',pr
plotimage,nhr,xst=4,yst=4,range=pr,title='HRDS HII Regions',min_dpi=100

ddr = mean(dr[1:*] - dr[0:*])
ddv = mean(dv[1:*] - dv[0:*])

cgAxis,xaxis=0,xr=[min(dr),max(dr)]+[-1,1]*ddr/2.,/xst,charsize=1.0,/save,$
       xtit='Physical Radius  [pc]'
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
cgAxis,yaxis=0,yr=[min(dv),max(dv)]+[-1,1]*ddv/2.,/yst,charsize=1.0,/save,$
       ytit='FWHM Velocity Dispersion  [km s!u-1!n]'
cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

cgContour,nhr,dr,dv,levels=findgen(21)*200.,/over,color=contcol

cgPlots,psym=14,color='brown',symsize=2,  16.5, 5.1
cgPlots,psym=4,color='black',symsize=2,   16.5, 5.1
cgPlots,psym=14,color='bisque',symsize=2,  7.4, 3.1
cgPlots,psym=4,color='black',symsize=2,    7.4, 3.1
cgPlots,psym=14,color='blu3',symsize=2,   14.7, 4.75
cgPlots,psym=4,color='black',symsize=2,   14.7, 4.75
cgPlots,psym=14,color='pink',symsize=2,   10.5, 3.6
cgPlots,psym=4,color='black',symsize=2,   10.5, 3.6
cgPlots,psym=14,color='grn4',symsize=2,    9.0, 4.1
cgPlots,psym=4,color='black',symsize=2,    9.0, 4.1

multiplot,/doxaxis,/doyaxis

;;================================
;; NEPXPO / NEPX

;; val = double(nepxpo) / double(nepx)
;; ;;val[where(~finite(val))] = 0.

;; print,m4_stat(nepxpo),m4_stat(nepx)
;; cgLoadct,13
;; pr = set_plot_range(val)
;; print,'NEPXPO / NEPX: ',pr
;; plotimage,val,xst=4,yst=4,range=pr,$;min_dpi=100,$
;;           tit='Parallax '+cgSymbol_extra('intersect')+' EMAF'

;; ddr = mean(dr[1:*] - dr[0:*])
;; ddv = mean(dv[1:*] - dv[0:*])

;; cgAxis,xaxis=0,xr=[min(dr),max(dr)]+[-1,1]*ddr/2.,/xst,charsize=1.0,/save,$
;;        xtit='Physical Radius  [pc]'
;; cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
;; cgAxis,yaxis=0,yr=[min(dv),max(dv)]+[-1,1]*ddv/2.,/yst,charsize=1.0,/save,$
;;        ytit='FWHM Velocity Dispersion  [km s!u-1!n]'
;; cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

;; cgContour,nepx,dr,dv,levels=findgen(21)*20,/over,color=contcol

;; cgPlots,psym=14,color='brown',symsize=2,  16.5, 5.1
;; cgPlots,psym=4,color='black',symsize=2,   16.5, 5.1
;; cgPlots,psym=14,color='bisque',symsize=2,  7.4, 3.1
;; cgPlots,psym=4,color='black',symsize=2,    7.4, 3.1
;; cgPlots,psym=14,color='blu3',symsize=2,   14.7, 4.75
;; cgPlots,psym=4,color='black',symsize=2,   14.7, 4.75

;; cgColorbar,range=pr,position=[!x.window[0],0.04,!x.window[1],0.05],$
;;            title='Well-Constrained %',charsize=0.8

multiplot,/doxaxis,/doyaxis

;;================================
;; F match


print,m4_stat(fma)
cgLoadct,13
pr = set_plot_range(fma)
print,'Matching: ',pr
plotimage,fma,xst=4,yst=4,range=pr,$ ;min_dpi=100,$
          tit='HRDS vs. EMAF + H!d2!n'

ddr = mean(dr[1:*] - dr[0:*])
ddv = mean(dv[1:*] - dv[0:*])

cgAxis,xaxis=0,xr=[min(dr),max(dr)]+[-1,1]*ddr/2.,/xst,charsize=1.0,/save,$
       xtit='Physical Radius  [pc]'
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
cgAxis,yaxis=0,yr=[min(dv),max(dv)]+[-1,1]*ddv/2.,/yst,charsize=1.0,/save,$
       ytit='FWHM Velocity Dispersion  [km s!u-1!n]'
cgAxis,yaxis=1,/yst,ytickformat='blank_axis'


cgContour,nma,dr,dv,levels=findgen(21)*50,/over,color=contcol

cgPlots,psym=14,color='brown',symsize=2,  16.5, 5.1
cgPlots,psym=4,color='black',symsize=2,   16.5, 5.1
cgPlots,psym=14,color='bisque',symsize=2,  7.4, 3.1
cgPlots,psym=4,color='black',symsize=2,    7.4, 3.1
cgPlots,psym=14,color='blu3',symsize=2,   14.7, 4.75
cgPlots,psym=4,color='black',symsize=2,   14.7, 4.75
cgPlots,psym=14,color='pink',symsize=2,   10.5, 3.6
cgPlots,psym=4,color='black',symsize=2,   10.5, 3.6
cgPlots,psym=14,color='grn4',symsize=2,    9.0, 4.1
cgPlots,psym=4,color='black',symsize=2,    9.0, 4.1

cgColorbar,range=pr,position=[!x.window[0],0.04,!x.window[1],0.05],$
           title='Distance Agreement %',charsize=0.8

print,!x.crange,!y.crange

myps,/done,/mp


END
